library(tidyverse)
library(tidybulk)
library(ggpubr)
library(ggrepel)

kls_count <- read_csv("Archive/FPP_TRPV3/data/Thanert-Pieper_2019_NSTI.csv", name_repair = make.names)

kls_meta <- read_delim("Archive/FPP_TRPV3/data/Thanert_meta.txt", name_repair = make.names) |>
  mutate(
    Sample = make.names(Sample),
    Location = replace_na(Location, "NA") |> make.names()
  ) |>
  select(c(Sample, Location, Classification))

kls_tibble <- kls_count |>
  pivot_longer(cols = where(is.numeric),
               names_to = 'Sample',
               values_to = 'abundance',
               values_transform = round) |>
  left_join(kls_meta)

kls_bulkobj <- kls_tibble |>
  filter(!is.na(Official.gene.symbol)) |>
  tidybulk(Sample, Official.gene.symbol, abundance)

kls_scaled <- kls_bulkobj |>
  aggregate_duplicates() |> # useful to merge transcripts in single gene, also can apply to data frame
  identify_abundant(factor_of_interest = Classification) |>
  scale_abundance()

# Plot count densities
kls_scaled %>%
  ggplot(aes(abundance_scaled, group=Sample, color=Classification)) +
  geom_density() +
  scale_x_log10()

# Reduce data dimensionality with arbitrary number of dimensions
kls_mds <- kls_scaled %>% reduce_dimensions(method="PCA", .dims = 3)

# Plot all-vs-all MDS dimensions 
kls_mds %>%
  pivot_sample() %>%
  ggplot(aes(PC1, PC2, colour=Classification)) +
  geom_point()

kls_adj <- kls_scaled |>
  filter(Classification %in% c('Others','Multiinfection')) |>
  adjust_abundance(~ Classification + Location) |>
  test_differential_abundance(~ Classification + Location)

# MA plot
kls_adj %>%
  keep_abundant() %>%
  pivot_transcript() %>%
  # Subset data
  mutate(significant = FDR<0.05 & abs(logFC) >=2) %>%
  mutate(feature = ifelse(significant, as.character(Official.gene.symbol), NA)) %>%
  # Plot
  ggplot(aes(x = logCPM, y = logFC, label=feature)) +
  geom_point(aes(color = significant, size = significant, alpha=significant)) +
  geom_text_repel() +
  scale_color_manual(values=c("black", "red")) +
  scale_size_discrete(range = c(0, 2)) 

kls_sig <- kls_adj %>%
  keep_abundant() %>%
  pivot_transcript() |>
  filter(FDR < 0.1)

# HIMR ------------
himr <- readxl::read_excel('mission/core_table_gene-HIMRSA-PBS.xls')
himr |> head() |> DT::datatable()

cytokine_of_interest <- c('Il6','Ccl20','Flt3l','Tslp','Csf2','Tnf','Tgfb','Il1b','Il2','Il4','Il10','Il17a','Il17f','Cxcl2','Cxcl10','Ccl2','Ccl4','Ccl17','Csf1','Csf3','Ifna4','Ifnb1','Infg')

gene_of_mva <- c('Hmgcr','Hmgcs1','Hmgcs2','Fdps','Mvd','Idi1','Idi2')

kegg_mva <- c(gene_of_mva, 'Pmvk','Acat1','Acat2','Ggps','Mvk')

gene_of_interest <- c(cytokine_of_interest, gene_of_mva)

tdb_himr <- himr |>
  select(gene_symbol...3, matches('^read_count')) |>
  pivot_longer(-1, names_to = 'sample', names_prefix = 'read_count_',
               values_to = 'count', values_transform = round) |>
  mutate(group = str_remove(sample, '\\d$')) |>
  rename(gene_symbol = gene_symbol...3)

tdb_himr |> write_csv('mission/FPP/himr_bulk.csv')

tdb_himr <- read_csv('mission/FPP/himr_bulk.csv')

tdb_himr <- tdb_himr |>
  tidybulk(.sample = sample, .transcript = gene_symbol, .abundance = count) |>
  aggregate_duplicates() |>
  identify_abundant(factor_of_interest = group,
                    .sample = sample,
                    .transcript = gene_symbol,
                    .abundance = count) |>
  scale_abundance()

# Plot count densities
tdb_himr %>%
  ggplot(aes(count_scaled, group = sample, color = group)) +
  geom_density() +
  scale_x_log10()

# Reduce data dimensionality with arbitrary number of dimensions
## MDS plot -------
tdb_himr %>%
  reduce_dimensions(method="MDS", .dims = 3) |>
  pivot_sample() %>%
  ggplot(aes(Dim1,Dim2,color = group, label = sample)) +
  geom_point() + geom_text_repel()

## pca plot ---------
tdb_himr |>
  reduce_dimensions(method="PCA", .dims = 2) |>
  pivot_sample() |>
  ggplot(aes(PC1,PC2,color = group, label = sample)) +
  geom_point() + geom_text_repel()

# use edgeR QLF by default
himr8h_adj <- tdb_himr |>
  filter(group != 'HIMR24hr') |>
  mutate(group = fct_relevel(group, 'PBS')) |> # relevel ref group
  test_differential_abundance(~ group)

himr24h_adj <- tdb_himr |>
  filter(group != 'HIMR8hr') |>
  mutate(group = fct_relevel(group, 'PBS')) |> # relevel ref group
  test_differential_abundance(~ group)

## volcano plot ----------
### 8h -----------
himr8h_sig <- himr8h_adj %>%
  keep_abundant() %>%
  pivot_transcript() 

himr8h_sig <- himr8h_sig |>
  mutate(type = case_when(FDR < .05 & logFC > 2 ~ 'Upregulated',
                          FDR < .05 & logFC < -2 ~ 'Downregulated',
                          .default = 'NS')) |>
  filter(!str_detect(gene_symbol, '^Gm'))

himr8h_named <- himr8h_sig |>
  filter(type != 'NS') |>
  group_by(type) |>
  slice_max(abs(logFC), n = 10)

# get description for gene symbol
himr8h_named |>
  select(gene_symbol) |>
  describe_transcript(.transcript = gene_symbol)

## highlight 10 highest logFC   
himr8h_sig |>
  ggplot(aes(logFC, -log10(FDR), color = type)) +
  geom_point() +
  scale_color_manual(values = c('deepskyblue','grey','orange')) +
  geom_vline(xintercept = c(-2,2), linetype = 'dashed') +
  geom_hline(yintercept = 1.3, linetype = 'dashed') +
  geom_text_repel(data = himr8h_named, aes(logFC, -log10(FDR), label = gene_symbol),
                           inherit.aes = FALSE, size = 5) +
  expand_limits(x = c(8,-8)) +
  theme_pubr() + labs_pubr() +
  ggtitle('8h after infection')

ggsave('mice_bulk_8h_volcano.pdf', width = 7, height = 6)

### FIG: 8h volcano ==========
## highlight panel of interest
himr8h_panel <- himr8h_sig |>
  filter(gene_symbol %in% c(key_cytokine, kegg_mva))

himr8h_sig |>
  ggplot(aes(logFC, -log10(FDR), color = type)) +
  geom_point(size = .3, alpha = .3) +
  scale_color_manual(values = c('blue','grey','red')) +
  geom_vline(xintercept = c(-2,2), linetype = 'dashed') +
  geom_hline(yintercept = 1.3, linetype = 'dashed') +
  geom_point(data = himr8h_panel, color = 'orange', size = 1) +
  geom_text_repel(data = himr8h_panel, aes(logFC, -log10(FDR), label = gene_symbol),
                           inherit.aes = FALSE, size = 1.5) +
  expand_limits(x = c(8,-8)) +
  theme_classic(base_size = 6, base_line_size = 1) +
  theme(legend.position = 'top', plot.margin = margin(), legend.box.margin = margin())

publish_pdf('micebulk/micebulk_8h_volcano.pdf')

### 24h -------
himr24h_sig <- himr24h_adj %>%
  keep_abundant() %>%
  pivot_transcript() 

himr24h_sig <- himr24h_sig |>
  mutate(type = case_when(FDR < .05 & logFC > 2 ~ 'Upregulated',
                          FDR < .05 & logFC < -2 ~ 'Downregulated',
                          .default = 'NS')) |>
  filter(!str_detect(gene_symbol, '^Gm'))

himr24h_named <- himr24h_sig |>
  filter(type != 'NS') |>
  group_by(type) |>
  slice_max(abs(logFC), n = 10)

## highlight 10 highest logFC   
himr24h_sig |>
  ggplot(aes(logFC, -log10(FDR), color = type)) +
  geom_point() +
  scale_color_manual(values = c('deepskyblue','grey','orange')) +
  geom_vline(xintercept = c(-2,2), linetype = 'dashed') +
  geom_hline(yintercept = 1.3, linetype = 'dashed') +
  geom_text_repel(data = himr24h_named, aes(logFC, -log10(FDR), label = gene_symbol),
                           inherit.aes = FALSE, size = 5) +
  expand_limits(x = c(11,-11)) +
  theme_pubr() + labs_pubr() +
  ggtitle('24h after infection')

ggsave('mice_bulk_24h_volcano.pdf', width = 7, height = 6)

### FIG: 24h volcano ===============
## highlight panel of interest
himr24h_panel <- himr24h_sig |>
  filter(gene_symbol %in% c(key_cytokine, kegg_mva))

himr24h_sig |>
  ggplot(aes(logFC, -log10(FDR), color = type)) +
  geom_point(size = .3, alpha = .3) +
  scale_color_manual(values = c('blue','grey','red')) +
  geom_vline(xintercept = c(-2,2), linetype = 'dashed') +
  geom_hline(yintercept = 1.3, linetype = 'dashed') +
  geom_point(data = himr8h_panel, color = 'orange', size = 1) +
  geom_text_repel(data = himr8h_panel, aes(logFC, -log10(FDR), label = gene_symbol),
                           inherit.aes = FALSE, size = 1.5) +
  expand_limits(x = c(8,-8)) +
  theme_classic(base_size = 6, base_line_size = 1) +
  theme(legend.position = 'top', plot.margin = margin(), legend.box.margin = margin())

publish_pdf('micebulk/micebulk_24h_volcano.pdf')

### FIG: mva dotplot =========
himr_panel <- himr24h_panel |>
  bind_rows(himr8h_panel,.id = 'test') |>
  mutate(test = case_match(test, '1' ~ '24 hpi', '2' ~ '8 hpi'))

himr_panel |>
  filter(gene_symbol %in% kegg_mva) |>
  ggplot(aes(test, gene_symbol,
             color = logFC, size = -log10(FDR))) +
  geom_point() +
  labs(x = 'Sample', y = 'Gene') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub

publish_pdf('micebulk/micebulk_mva_bubble.pdf')

### FIG: cytokine dotplot ========
himr_panel |>
  filter(gene_symbol %in% key_cytokine & str_starts(test, '8')) |>
  ggplot(aes(test, gene_symbol,
             color = logFC, size = -log10(FDR))) +
  geom_point() +
  labs(x = 'Sample', y = 'Gene') +
  scale_color_gradient2(low = 'blue', high = 'red') +
  theme_jpub

publish_pdf('micebulk/micebulk_8h_cytokine_bubble.pdf', width = 40)

## save result data ----
himr24h_sig <- himr24h_adj %>%
  keep_abundant() %>%
  pivot_transcript()

himr24h_sig |>
  bind_rows(himr8h_sig,.id = 'test') |>
  mutate(test = case_match(test, '1' ~ '24h vs PBS', '2' ~ '8h vs PBS')) |>
  select(gene_symbol, logFC, FDR, test) |>
  write_csv('mission/himrsa-pbs-logFC-edgeR-QLF.csv')

tdb_himr |>
  select(gene_symbol, sample, group, count_scaled) |>
  write_csv('mission/himrsa-pbs-TMM-scaled_count.csv')

himr24h_sig <-
  read_csv('mission/himrsa-pbs-logFC-edgeR-QLF.csv') |>
  filter(str_starts(test, '2'))

himr8h_sig <-
  read_csv('mission/himrsa-pbs-logFC-edgeR-QLF.csv') |>
  filter(str_starts(test, '8'))
